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Abstract. We present and analyze a micro/macro acceleration technique 
for the Monte Carlo simulation of stochastic differential equations (SDEs) in 
which there is a separation between the (fast) time-scale on which individual 
trajectories of the SDE need to be simulated and the (slow) time-scale on 
which we want to observe the (macroscopic) function of interest. The method 
performs short bursts of microscopic simulation using an ensemble of SDE 
realizations, after which the ensemble is restricted to a number of macroscopic 
state variables. The resulting macroscopic state is then extrapolated forward 
in time and the ensemble is projected onto the extrapolated macroscopic state. 
We provide a first analysis of its convergence in terms of extrapolation time 
step and number of macroscopic state variables. The effects of the different 
approximations on the resulting error are illustrated via numerical experiments. 



1. Introduction 

In many applications, one considers a process that is modeled as a stochastic 
differential equation (SDE), while one is ultimately interested in the time evolution 
of the expectation of a certain function of the state, i. e., in weak approximation. 
Consider, for instance, the micro/macro simulation of dilute solutions of polymers 



34 , which will be the motivating example in this paper. Here, an SDE models 
the evolution of the configuration of an individual polymer driven by the flow field, 
and the function of interest is a non-Newtonian stress tensor (the expectation of a 
function of the polymer configuration) . For this type of problem, one often resorts 
to Monte Carlo simulation |6j, i. e., the simulation of a large ensemble of realizations 
of the SDE, combined with ensemble averaging to obtain an approximation of the 
quantity of interest at the desired moments in time. For concreteness, we introduce 
the SDE 

(1) dX(t) = a(t,X(t)) dt + b(t,X{t))*dW(t), t e I := [t° ,T], X(t°) = X Q , 

in which a : I x R d -> R d is the drift, b : I x R d -> R dxm is the diffusion, and 
W(t) is an m-dimcnsional Wiener process. The initial value Xq is independent of 
W and follows some known distribution with density tpo(x). As usual, (JT|) is an 
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abbreviation of the integral form 

X(t) = X Q + f a(s,X(s)) ds + f b(s,X(s))*dW{s), tel. 
Jt° Jt° 

The integral with respect to W can be interpreted, e.g., as an Ito integral with 
* dW(s) — dW(s) or as a Stratonovich integral with * dW(s) = o dW(s). The 
function of interest for the Monte Carlo simulation is defined as the expectation E 
of a function f(X(t)j, 

f(t)=Ef(X(t)). 

The numerical properties of Monte Carlo simulations have been analyzed ex- 
tensively in the literature. We mention studies on the order of weak convergence 
of explicit |7)[l2)[30)[3lJ[40)[48] and implicit j8pl|ll||30H32] time discretizations of 



SDE ([lj, the investigation of stability 19 20 22 49 56 , and techniques for variance 
reduction |16[|17[|44| . For more references, we refer to 30 . Also for strong approxi- 
mation, there has been a growing interest in the study of numerical methods for 
stiff SDEs [T^[5}[30j[4T}[42j[^}[55] . 

In this paper, we present and analyze a micro/macro acceleration technique 
for the Monte Carlo simulation of SDEs of the type in which there exists a 
time-scale separation between the (fast) time-scale on which individual trajectories 
of the SDE need to be simulated and the (slow) time-scale on which the function /(f) 
evolves. The proposed method is motivated by the development of recent generic 
multiscale techniques, such as equation-free [ 28 , 29 and heterogeneous multiscale 
methods [I3||l4] . We use the simulation of a dilute polymer solution as an illustrative 
example. The microscopic level is defined via an ensemble X = (Xj)j =1 of J 
realizations evolving according to equation ([I]); the macroscopic level will be defined 
by a set of L macroscopic state variables U = (f/jOjLii with Ui(t) = FiUi(X(t)), 
for some appropriately chosen functions ui. The method exploits a separation in 
time scales by combining short bursts of microscopic simulation with the SDE 
([lj with a macroscopic extrapolation step, in which only the macroscopic state U 
is extrapolated forward in time. One time step of the algorithm can be written 
as follows: (1) microscopic simulation of the ensemble using the SDE 0; (2) 
restriction, i. e., extraction of (an estimate of) the macroscopic state (or macroscopic 
time derivative); (3) forward in time extrapolation of the macroscopic state; and (4) 
matching of the ensemble that was available at the end of the microscopic simulation 
with the extrapolated macroscopic state. Remark that the resulting method is fully 
explicit as soon as the microscopic simulation is explicit, and that the method can 
readily be implemented as a higher order method by an appropriate choice of the 
extrapolation. 

The main contributions of the present paper are the following: 

• From a numerical analysis viewpoint, we study convergence of the proposed 
micro/macro acceleration method in the absence of statistical error. Specif- 
ically, we discuss how the error that is introduced during the matching 
depends on the number of macroscopic state variables, and how the deter- 
ministic error depends on the extrapolation time step. Additionally, we also 
comment on the effects of the extrapolation on the statistical error. Finally, 
we give a general convergence result. 

• From a practical viewpoint, we provide numerical results for a nontrivial 
test case, showing the interplay between the different sources of numerical 
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error. We illustrate the effects of the choice of macroscopic state variables, 
as well as the dependence of the numerical error on the chosen extrapolation 
strategy. 

A stability analysis of the proposed method will be given in a separate publication. 

The remainder of the paper is organized as follows. In Section[2] we first introduce 
the mathematical setting of the paper. We discuss the necessary assumptions on the 
SDE for our Monte Carlo setting, and introduce the illustrative example that 
will be used for the numerical experiments. In Section |3j we propose the algorithm 
that will be the focus of this paper. Section [4] provides some results on the matching 
operator, whereas the extrapolation operator is discussed in Section [5j We provide 
a general convergence result in Section [6] Section [7] provides numerical illustrations, 
which are chosen to illuminate the properties of the proposed method. We conclude 
in Section [8] where we also outline some directions for future research. 



2. MATHEMATICAL SETTING 

In this section, we introduce in detail the notations that we will use (Subsec- 



tion 



2.1 1, as well as the illustrative example that will be considered throughout the 



paper (Subsection 2.2 1 



2.1. Notations. We first define the appropriate function spaces. Let Cp(E d ,M) 
denote the space of all g £ C r (M. d ,R) fulfilling that there are constants C > and 
k > such that \dl.g(x)\ < (7(1 + ||cc||' t ) for any partial derivative of order i < r and 
all x € R d . Further, let g € C Q p r (I x R d ,R) if g(-,x) € C q (I,R), g{t, •) € C7 r (K d ,K) 
for all t £ / and x £ R d , and |<9^g(s, a:)| < C(l + |Ur|| K ) holds for < i < r uniformly 
with respect to s S [t°, t] and for all x £ R d |oj 40 . 

We consider the SDE 0. Besides the exact solution X(t) of SDE starting 
from X(t°) = Xq, we also introduce the exact solution of an auxiliary initial value 
problem for the SDE: the solution of the SDE starting from an initial value Z at 
time t will be denoted as X*' Z . With this notation, we state the following definition: 

Definition 2.1 (Uniform weak continuity of the SDE). Consider a class of random 
variables. The SDE is called uniformly weakly continuous for this class if, for all 
g £ Cp 2 (I x there exist constants Ato > and C such that 

(2) \Eg(s,X t ' z (t + At))-Eg(s,Z)\<CAt 

holds for all initial values Z in the considered class, At £ [0,Aio], and all t £ 
[t°,T - At], s € I. 

Next, we discretize equation ([T]) in time with step size St and denote the numerical 
approximation Y k — Y(t k ) « X(t k ) with t = t° + kSt. Leaving the concrete 
choice of discretization undecided for now, we introduce a short-hand notation for 
an abstract one step discretization scheme, 

(3) Y k+1 = 4>{t k , Y k ;St), Y° = X , k> 0. 

In the Monte Carlo setting, we are interested in weak approximation of the SDE 
0. We define the weak order of consistency of the discretization <f) as follows 
(compare |40]): 



4 



KRISTIAN DEBRABANT AND GIOVANNI SAMAEY 



Definition 2.2 (Weak order of consistency of SDE discretization). Assume that 
for all g G Cp' 2(p * +13 (/ x R d ,R) there exists a C g € C P (R d ,E) such that 

\Eg(s,i>(t,Z;6t)) -Eg(s,X t ' z (t + 5t))\ < C g (Z) St"* +1 

is valid for Z S ]R d , s G J, and t, t + St e [t°,T\, Then the one step method is 
called weakly consistent of order . 

Since we can only simulate a finite number of realizations of approximations of 
([lj (via its discretization ([3])), we also need to approximate the expectation E by 
an empirical mean E. Using the ensemble y = (Yj)j =1 of realizations of (|3]), the 
empirical mean E is defined as 

E/(y) = ~x;/(is)- 

J 3 = 1 

The numerical integration scheme for the ensemble y will be denoted as 

(4) y k +i = (h (t k ,y k ;5t)= (cj>{t\y);8t)) J ^. 

The total error of a Monte Carlo simulation of (|T|) consists of a deterministic 
error due to the time discretization, and a statistical error due to the finite number 
of realizations. In general, for a given tolerance, a Monte Carlo simulation will be 
most efficient if the statistical and deterministic error are balanced. However, as 
for stiff systems of ODEs, one might encounter situations in which the time step St 
cannot be increased because of stability problems. The acceleration method that 
we will propose in Section [3] will prove to be particularly useful to accelerate Monte 
Carlo simulations in such situations. 

Remark 2.3 (Probability density functions). One can equivalently describe the pro- 
cess via an advection-diffusion equation, also known as Fokker-Planck equation 
(see, e.g., (47]), which in the Ito case takes the form 

(5) e> t¥ >=-V(a^)+*V- [V • (b T b^)} , 

and which describes the evolution of the probability density function <p(t, x) of X(t), 
starting from the initial density <p(t , x) = (po(x). The functional of interest then 
becomes 

(6) /(t) = J f(xMt,x)dx. 

For simulation purposes, however, the Monte Carlo algorithm is generally preferred, 
due to the possibly high number of dimensions of the Fokker-Planck equation. 

Remark 2.4 (Spatial dimension). While the model ([lj can have arbitrary dimension, 
we will consider a one-dimensional version in the numerical illustrations for ease 
of visualization. In that case the state vector reduces to a scalar, X(t). Whenever 
we consider the one-dimensional case, all bold typesetting in equation (JlJ will be 
removed. 
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2.2. A motivating model problem: FENE dumbbells. To illustrate the behav- 
ior of the proposed numerical methods, we will consider the micro/macro simulation 
of the evolution of immersed polymers in a solvent. Here, one models the evolution 
of the configuration of a polymer ensemble via an SDE of the type ([I]) , driven by 
the flow field, for each of the individual polymers. This results in a polymer configu- 
ration distribution at each spatial point. This microscopic model is coupled to a 
Navier-Stokes equation for the solvent, in which the effect of the immersed polymers 



is taken into account via a non-Newtonian stress tensor. We refer to [23 34 36 for 
an introduction to the literature on this subject. 

In this paper, we consider only the Monte Carlo simulation of the microscopic 
model, leaving the coupling with the Navier-Stokes equations for future work. We 
eliminate the spatial dependence by considering the microscopic equations along the 
characteristics of the flow field, i. e., in a Lagrangian frame. In general, a microscopic 
model describes an individual polymer as a series of beads, connected by nonlinear 
springs, resulting in a coupled system of SDEs for the position of each of the beads. 
In the simplest case, that we will also use as an illustrative example here, one 
represents the polymers as non-interacting dumbbells, connecting two beads by 
a spring that models intramolecular interaction. The state of the polymer chain 
is described by the end-to-end vector X(t) that connects both beads, and whose 
evolution is modelled using the non-dimcnsionaliscd SDE 



(7) dX(t) = 



dt+ ^=dW(t), 
VWe 



where n(t) is the velocity gradient of the solvent, We is the Weissenberg number, 
and F is an entropic force, here considered to be finitely extensible nonlinearly 
elastic (FENE), 

with 7 a non-dimensional parameter that is related to the maximal polymer length. 
The resulting non-Newtonian stress tensor is given by the Kramers' expression, 

(9) r p (t) = ^(E(x(t)»F(X{t)))-Idj, 



in which e represents the ratio of polymer and total viscosity; see [36] for details and 
further references. This model, which is of the type ([lj, takes into account Stokes 
drag (due to the solvent velocity field), intramolecular elastic forces, and Brownian 
motion (due to collisions with solvent molecules) . The functional of interest in the 
Monte Carlo simulation is f(t) — T p (t). 

Equations ^ and ^ ensure that the length of the end-to-end vector, \\X\ 



cannot exceed the maximal value y/j 26 . However, a naive explicit discretization 



scheme might yield polymer lengths beyond this maximal value. This can be avoided 
via an accept-reject strategy, see, e.g., [45] Section 4.3.2]. Here, for each polymer, 
the state after each time step is rejected if the calculated polymer length exceeds 



(1 — VSt)j, and a new random number is tried until acceptance. To prevent the 
distribution of the approximation process to be heavily influenced, the microscopic 
time-step has to be chosen small enough. Alternatively, one can use an implicit 



method 45 , which alleviates the time-step restriction. However, even for implicit 



SDE discretizations, the maximal time step is limited when coupling the Monte Carlo 
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simulation with a discretization of the Navier-Stokes equations for the solvent. This 
is due to the fact that the coupling between the Monte Carlo and Navier-Stokes parts 
is, in most existing work, done explicitly in time, creating an additional stability 
constraint due to the coupling. (Some notable exceptions are given in |35||52| .) With 
the algorithm that we propose, one would extrapolate both the Monte Carlo and 
the Navier-Stokes part of this coupled simulation simultaneously. This will be done 
in future work. Here, we simply conclude that, for this model problem, the required 
time step for a stable SDE (or coupled) simulation may indeed be small compared 
to the time scale of the evolution of the stress. 



As in, e.g., 27 , we will consider a one-dimensional version in the numerical 
illustrations; the stress tensor then reduces to a scalar r p (t). As time discretization 
we will use the explicit Euler-Maruyama scheme, combined with an accept-reject 
strategy. 

3. Micro /macro acceleration method 

In this section, we describe the micro/macro acceleration algorithm that is the 
focus of the present paper. As said above, the goal of the method is to be faster than 
a full microscopic simulation, while converging to a full microscopic simulation when 
the extrapolation time step vanishes (see Section [6]) . The algorithm combines short 
bursts of microscopic simulation of an ensemble X of J realizations of the SDE ([lj 
with a macroscopic extrapolation step. During this macroscopic extrapolation, only 
a set of L macroscopic state variables U are extrapolated forward in time, and the 
microscopic ensemble then needs to be matched onto the extrapolated macroscopic 
state. 

3.1. Description of algorithm. Before giving a detailed description of the algo- 
rithm, let us first elaborate slightly on the macroscopic state variables, U = {Ui)^ =1 , 
which arc defined as expectations of scalar functions ui of the state X and time t, 

(10) U l (t) = Eu l (t,X(t)). 

Remark 3.1 (Choice of macroscopic state variables). The choice of the functions 
ui is problem-dependent and will be specified with the numerical illustrations for 
the examples considered in this text. For the exposition in this section, it may be 
helpful to think about the standard moments of the distribution in a one-dimensional 
setting, i.e., ui(t,x) — x l . Note that by allowing u\ to depend directly on time t, 
centralized moments u\{t,x) = x, Ui(t,x) = {x — U\(t)) 1 for I > 2, for instance, can 
also be considered. 

We introduce two abstract operators to connect both levels of description: a 
restriction operator 

(11) n-.y^u = n(y), 

that maps a microscopic ensemble onto a macroscopic state, and a matching operator 

(12) v : u,y* h> y = v(u,y*), 

which matches a given microscopic ensemble y* with an imposed macroscopic state 
U. Both operators will be discussed in detail in Subsection |3.2| 

We further introduce some notation. Let I At — {t°, t 1 , . . . , t N } be a (macroscopic) 
time discretization of the time interval /, with t° < t 1 < . . . < t N = T and step 
sizes At n = t n+1 — t n for n = 0, 1, . . . , N — 1. We also introduce the discrete time 
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instances t n ' — t n + kSt that are denned on a microscopic grid, and, correspondingly, 
the discrete approximations U n,k w U(t n,k ) and U n w U(t n ), and analogously for 
y and Y. Clearly, (•)" = (-)"' . 

One step of the micro/macro acceleration method then reads: 

Algorithm 3.2 (Micro/macro acceleration). Given a microscopic state y n at time 
t n , advance to a microscopic state y n+1 at time t n+1 via a three-step procedure: 

(i) Simulate the microscopic system over K time steps of size St, 

y n ' k = 0y(t n '*~ 1 ,y n,fc ~ 1 ;&), l<k<K, 
and record the restrictions U n,k — lZ(y n ' k ), as well as an approximation of 
the function of interest f n > k = E/(y™' fc ). 

(ii) Extrapolate the macroscopic state U from (some of) the time points t l ' k ,i = 
0, . . . , n, k = 1, . . . , K, to a new macroscopic state U n+1 at time t n+ , 

(13) U"+^£((U- k ):^ 00 -(AU)- =0 ,St). 

(iii) Match the microscopic state y n,x with the extrapolated macroscopic state 

yn+1 _ -p(yn,K ^ Jjn+ly 

Some basic requirements for the matching and restriction operators are given 
in Subsection |3.2| Specific choices and convergence properties for the algorithmic 
components are given in Section [3] (matching) and Section [5] (extrapolation). 

Remark 3.3 (Closure approximation). Due to the possibly high computational 
cost of Monte Carlo simulation, another route has been followed in the literature, in 
which one derives an approximate macroscopic model to describe the system; see, 
e.g., 1 2 1 , 24 ,27,38,3951 for derivations of macroscopic closures for FENE dumbbell 
models. The goal then is to obtain a closed system of L evolution equations for the 
macroscopic state variables U , complemented with a constitutive equation, for the 



observable / of interest as a function of these macroscopic state variables. In 25 , a 
quasi- equilibrium approach is proposed, based on thermodynamical considerations; 
although the method has been formulated for the FENE dumbbell case, it is 
applicable to general SDEs of the type 0. Several algorithms have been presented 



to simulate the evolution of the quasi-equilibrium model numerically 50 59 . Note 
that, in contrast with the method presented in the present paper, numerical closures 
introduce the modeling assumption that a closed model in terms of the macroscopic 
state variables exists. The micro/macro acceleration method presented here only 
uses these macroscopic state variables for computational purposes, and maintains 
convergence to the full microscopic dynamics (see Section [6]). 

3.2. Matching and restriction operators. Next, we give some detail on the 
restriction and matching operators that connect the microscopic and macroscopic 
levels of description. To obtain the macroscopic state from the microscopic ensemble, 
the restriction operator is defined, which can readily be obtained by replacing the 
expectation E by the empirical mean E, 

(14) u l (t) = n l (y(t)) = Eu l (t,y(t)). 

Due to the constraint lZ(y n+1 ) = U n+1 , during the matching step, the elements 
of y n+1 are in general not independent, and Monte Carlo error estimates are not 
straightforward to obtain. For example, the Central Limit Theorem is not directly 
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applicable [6]. Therefore, to discuss the matching operator V, we first turn to the 
idealized operators V and 1Z obtained in the limit J — > oo, eliminating statistical 
error. Rather than acting on ensembles of configurations, these operators are defined 
directly on the random variables. More precisely, the restriction operator 1Z reduces 
a random variable Z to macroscopic state variables, 

(15) n(Z) = {Ki{Z)) L l=l with %{Z) = Ui =Eui(Z) for I = l,...,L. 

The matching operator V maps a random variable Z onto a new random variable 
Z* that corresponds to a macroscopic state U* , 

(16) Z* =V(Z,U*) with K(Z*) = U*. 

When analyzing the matching, we will consider the idealized operators V and 1Z. 

A priori, we allow for considerable freedom in the definition of the matching 
operator, requiring only a few properties to be satisfied for any reasonable matching. 
First, we impose a projection property: 

Definition 3.4 (Projection property). A matching operator V associated to a 
restriction operator 1Z satisfies the projection property if 

U = TZ(Z) => Z = V[Z, U) = Vu(Z) 

for any suitable random variable Z, i.e., if Vjj is a projection operator for every U. 

The projection property states that a random variable remains unaffected by 
projection if its macroscopic state is equal to the m acro scopic state on which one 



wants to project. We remark also that Definition 3.4 implies Vfj = Vjj, which 
justifies the use of the term ■projection. 

Next, we consider the number of macroscopic state variables L to vary, and 
define a sequence of vectors of macroscopic state variables \U[l]) L _ x 2 , such that 

U[l] = (Ui)f =1: i.e., for increasing L, additional macroscopic state variables are 
added. The corresponding sequences of matching and restriction operators are 
denoted as (V[l]) l=1 2 and (T^[l)) L=1 2 , respectively. 

Using this notation, we are ready to formulate the definitions of continuity and 
consistency of the matching step: 

Definition 3.5 (Continuity of matching). Consider a set of random variables 
and a set of sequences of macroscopic states. A sequence of matching operators 
{V[l]) l1 2 is called continuous for these sets if, for all g £ Cp 2 (I x K d ,K), there 
exists a constant C, depending only on g, such that 

(17) \Eg(s,V [L] (Z,U^)) -E 5 ( S) P m (Z,Z7+))| < C\\U* [L] - U+ L] \\ 

holds for all L > 1, all sequences (u* L -^J , (jJ^^j of macroscopic 

states and all Z in the considered sets, and all s E I. 

Here and in the following, the norm || • |j should be chosen such that it remains 
bounded for L — > oo for component-wise bounded sequences. 

Definition 3.6 (Consistency of matching). Consider a sequence of matching op- 
erators [P[l]) r_ x 2 ■ This sequence is called consistent for a class of sequences 

of triples (^Zy L y ^[L]j 01 random variables Z* L ^, anc ^ macroscopic 
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states (U[l]), if for all g £ Cp 2 (I x E d ,K) there exist constants Cl, with Cl -> 

for L — y oo, and I/o such that it holds 

(18) 

\Eg( S ,V [L] (Zfa,U [L] ))-Eg(s,P [L] (Z+ L] ,U [L] ))\ < C L \E g(s, Zfa) -E g(s, Zfa )| 

for all L > Lq and all sequences of triples {zfa, Z^, U[L]j in the considered 

class. 

The possible dependence of the random variables on L is present since the 
random variables will be considered to have been generated using the micro/macro 
acceleration algorithm, and therefore depend on L. 

Whereas continuity measures (in a weak sense) the difference between the match- 
ing of a random variable with two different macroscopic states, the consistency 
measures the difference between the matching of two different random variables 
with the same macroscopic state. 

Definitions |3.4[ |3.6| and |2.1| immediately imply the following corollary: 

Corollary 3.7. Consider a sequence of matching operators (^ 5 [L]) i _ 12 > an< ^ 
assume that this sequence is consistent for a set of sequences of triples 

(z m , X*~' z m (t*),K [L] (X*~> z m (t*))) l=i a , 

where Z\n are random variables, t~ = t* — At, At £ [0,t* — t ], t* £ /. Suppose 
further that SDE ([l]) is uniformly weakly continuous for the set of all Z[ L ] . Then 
for all g £ Cp 2 (I x R d ,R) there exist constants Ato > and Cl, with Cl for 
L — > oo, and Lq, such that it holds 

(19) \Eg(a,X t ~' z W{t*))-Eg( 8 ,V lL] (Z lL] ,U^)\ < C L At 

for all L > L , all (z [L] , X l ~ > z w (t* ), Ufa =K [L] (X t ~ < z m(t*)f) in the 

considered set, At £ [0, At ], and all t* £ [t° + At, T], s £ I. 

This corollary states that the difference, measured in a weak sense, between the 
exact distribution at some time instance t* and a matching from a previous time 
t~ = t* — At with the exact macroscopic state at t* vanishes when matching with 
more macroscopic state variables or letting At — > 0. 

In the remainder of the text, we will only use the properties stated above. 
Therefore, one may use any matching operator that would come to mind for a 
particular problem, as long as the above properties are satisfied. 

3.3. Matching failure and adaptive time-stepping. In the numerical experi- 
ments, one might encounter situations in which the distributions evolve on time-scales 
that are not significantly slower than those of the macroscopic functions of interest. 
In that case, when taking a large extrapolation time step, the extrapolated macro- 
scopic state differs significantly from the state corresponding to the last available 
polymer ensemble, and it is possible that matching the ensemble with the desired 
macroscopic state fails. Consequently, a failure in the matching can be used as an 
indication to decrease the extrapolation time step. Based on this observation, we 
propose the following criterion to adaptively determine the macroscopic step size 
At: If the matching fails, we reject the step and try again with a time step 

(20) At ncw = max(aAt,KSt), a < 1, 
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whereas, when the matching succeeds, we accept the step and propose 
(21) Ai new = min(aAi, Ai max ), a > 1 

for the next step. If the macroscopic step size A< now = KSt, matching becomes 
trivial (the identity operator), since there is no extrapolation. Note that, when this 
happens, the criterion will ensure that the larger time steps are tried after the next 
burst of microscopic simulation. 

4. Matching operator 
In this section, we analyze the matching operator in more detail. We first define 



the specific matching operator that will be considered in this paper (Subsection 4.1 ) 



We then prove a consistency result in a very special case (Subsection 4.2 1 and give 
and illustrate then a conjecture for the general case (Subsection |4.3[). 



4.1. Definition of specific matching operator. The matching operator V can 
be defined in several ways. One option is to minimize the difference between the 
given ensemble and the result of the matching, for instance in the 2 — norm, i. e., 

(22) y n+1 = argmin hz-y n > K \\ 2 . 

Z-. n(z)=u n ^ * 

An approximation of the resulting ensemble can be obtained using a Lagrange 
multiplier technique: 



(23) 



L 

yn+l = yn,K + XtVyUtiy 71 '*), 
1 = 1 



with A = {A;}f =1 such that Ki{y n+1 ) = Up +1 for I = 1, . . . , L. 

Then, y n+1 — V(y n ' K , U n+1 ) results after obtaining A from a Newton procedure 
that solves the L-dimensional nonlinear system that defines the constraints. (Note 



that one can show that the resulting ensemble exactly satisfies (22 1 if we write an 
implicitly defined gradient Vy7Zi(y n+1 ) in the first line of (231.) 



By applying the implicit function theorem, we obtain the following lemma: 



Lemma 4.1. The problem (23) has a locally unique solution if 

• in the case of standard empirical moments Up , 
(24) det([/? +fc _ 2 ). fc=l! L ^0, 

• in the case of empirical centralized moments Up , 

Remark 4.2. The determinant 

det(EJ^ a (t»)) 



obtained by replacing the empirical moments in (24 1 by the moments themselves, 
is, from the theory of the Hamburger moment problem, known to be positive for 
distributions with finite support. 

Alternative matching operators can be defined by measuring the difference be- 
tween random variables differently, for instance using the Kullback-Leibler divergence 
(relative entropy) [33] . For some choices of the macroscopic state variables U one 
could also consider classical moment matching, as described in [6]. 
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Remark 4.3 (Matching for the FENE dumbbells). For the FENE dumbbells, an 
accept-reject strategy is applied during the combined evolution and matching, i.e., 
if the state of a polymer would become unphysical during the matching, we reject 
the trial move for that specific polymer in the evolution step and repeat the time 
step for this polymer, after which the matching of the ensemble is tried again. 

4.2. Particular results for normal distributions. For simplicity of the argu- 
ment, and without loss of generality, we restrict ourselves to a one-dimensional 
notation for the rest of this section. In the case of scalar normally distributed 
random variables, the following result can be proved. 

Lemma 4.4. Consider a scalar random variable Z and suppose that the macroscopic 
state variables are given by U\(z) = z and u<2,{z) — (z — Ui) 2 . Suppose further that 



in analogy to (22 \, V is given by 



(25) V(Z, U*) = argmin ]- E ((Z* - Z) 2 ) . 

Then, the random variable V(Z,U*) is also normally distributed with mean and 
variance given by U* . 



Proof. Denote by 



T 2 



Then, (251 yields 



(26) V(Z,U*)=±JK^\Z-n±»* i 



Thus, if Z is normally distributed, so is V(Z, U*). □ 

With the help of this lemma, we can easily show the following two corollaries. 

Corollary 4.5. Suppose that the hierarchy of macroscopic state variables is defined 
using Ui{z) — z and ui(z) = {z — Ui) 1 for I > 2. Suppose further that V is given 



by (251. Then, the sequence of matching operators is consistent with Cl = for 
L > 2 = L for all sequences of triples {z^ L y ^[L\i ^[L]j > w ^ ere an< ^ ^[L] 

are normally distributed and (u^ L ^j are sequences of centralized moment 

values consistent with normal distributions. 

Proof. We first consider the case L — 2. As the normal distribution is uniquely 



determined by its first two (centralized) moments, Lemma 4.4 implies that for two 
normally distributed random variables Z\ and Zi, V\^\(Z\, UrL) and 7 3 j 2 ]('^2, ^[2]) 



are identically distributed, and thus C2 = 0. For the same reason, (26 1 holds also 



for L > 2, and also in this case Cl = 0, and the matching is consistent. □ 

Corollary 4.6. Suppose that the hierarchy of macroscopic state variables is defined 
using u\{z) — z and ui(z) = (z — U\) 1 for I > 2. Suppose further that V is given 



by (25 I. Then, the sequence of matching operators is continuous for all normally 
distributed random variables, and sequences of centralized moment values consistent 
with normal distributions. 



12 



KRISTIAN DEBRABANT AND GIOVANNI SAMAEY 



Proof. For L > 2, Lemma 4.4 implies again that Z* = V\l]{Z,U? l ^) and Z + = 
V[l](Z, U^t) are normally distributed if Z is normally distributed and U^ L y t/r^, 
are sequences of centralized moment values consistent with normal distributions. 
Denoting the corresponding expectations by /i*, resp. [i + , and variances by (c*) 2 , 
resp. (cr+) 2 , it holds for all / <= C£(Kt,R) 



|E/(*,Z*)-E/(s,Z+)| = 



(/(«, (J*z + (j,*) - f(s, a+z + -Le" z2/2 dz 



2tt 



V Z7T 



where £ z <G [a*z + fj,*, a + z + As there exist constants C and r £ N such that 
Cz)l < ^(1 + |max{cr*,cr + }z + max{/j*, ^ + }| 2r ), this implies also that the 



sequence of matching operators is continuous (the corresponding equation (17 1 can 



be verified similarly in the case L = 1). □ 

Several observations can be made. First, one can obtain a similar result whenever 
the distributions are defined by a finite number of moments (for instance, a lognormal 
distribution) by taking this knowledge into account when defining the matching 
operator. Second, we remark that, if the random variables are normally distributed 
at all moments in time, this implies that they represent solutions of a linear SDE in 
the narrow sense, 

(27) dX{t) = (ai(t)X(t) + a 2 {t)) dt + b(t) * dW{t), t € I, 

with normally distributed initial values. For this equation, it is clear that the 
complete time evolution of the distributions can be completely described by a 
system of two ODEs for U\ = \x = E X and U% = a 1 = E ( (X — /i) 2 ) , namely 



(28) 



dUi/dt = ai(t) Ui +o 2 (t), 
dU 2 /dt = 2a x {t) U 2 + b(t) 2 . 



As a consequence, the matching operator (25 1 corresponds to the reconstruction of 



the normal distribution corresponding to the given macroscopic state. 

4.3. Conjecture for general distributions. We now turn to more general dis- 
tributions. We assume that the distribution of the random variable is uniquely 
determined by its moments; this is the case if, for instance, the moment generating 
function V" 00 ,, E ^, ^ is bounded in an interval around 0. 
In this setting, we propose the following conjecture: 

Conjecture 4.7. Consider a sequence of restriction operators (^■[i]) i _ 12 m 
which the macroscopic state variables corresponding to 7Z[l] & r e defined as the 
first L centralized moments of the distribution, i.e., u\(z) — z, Ui(z) = (z — Ui) , 
for / = 2,...,i, and define the corresponding sequence of matching operators 
vl L ]) L-i 2 v ^ a '25). Consider further a set S of random variables for which all 



(centralized) moments of its members exist and uniquely determine the corresponding 
distribution function, and each moment can be uniformly bounded. Then, the 
sequence of matching operators is continuous for S and all sequences of macroscopic 
states Ujx], and consistent for all sequences of triples {z* L -^ , Z^ , £/[z,]J where 

Z [LY Z [L\ E S - 
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(a) Empirical probability density functions (b) Error of the i-th moment as a function of I 



FIGURE 1. Results after projecting a prior ensemble of FENE 
dumbbells onto the first L even centralized moments of a reference 
ensemble for several values of L. Simulation details are given in 
the text. 



We illustrate the main properties of the matching operator for general distributions 
by means of numerical experiments. Below, we consider equation Q in one space 
dimension, with n{t) = 2, F(X) the FENE force (|8) with 7 = 49, We = 1. 
We discretize in time with the classical Euler-Maruyama scheme with time step 
8t = 2 • 10~ 4 , and simulate J = 1 • 10 5 realizations, whose initial state at time 
t° = is taken from the invariant distribution of equation ^ for n(t) = 0. As the 
macroscopic state U[lU we consider the first L even centralized moments, since for 
the exact solution, the odd (centralized) moments vanish due to symmetry. 

4.3.1. Error dependence on the number of moments. We first simulate up to time 
t* = 1.15 and record the microscopic state y* and corresponding macroscopic states 
U* L ^ for L = 1, . . . , 10, at time t* , as well as the microscopic state y~ at time 

t~ = 1. We then project the ensemble y~ onto the macroscopic state t/j^j and 

compare the density of "P[L](y - , w ith that of y* . Since the absolute value 
of the moments increases quickly with the order of the moment, the residuals in 
the Newton procedure for the matching are scaled relative to the requested value of 
the corresponding moment; the Newton iterations are stopped if the norm of the 
residual is smaller than 1 • 10 -9 . 

We perform three tests to examine the convergence (in empirical density) of 
V[L](y~ , r^jx]) t° y* f° r L — > 00. First, we visually inspect the corresponding 
empirical probability density functions, see Fig. |l(a) Shown are histogram 



approximations of the empirical density of \y~\ (the initial condition for the 
matching), the reference empirical density (p* of \y*\, and approximations (p^ 
of V[L](y~, Um) for several values of L. The figure visually suggests that, when 
increasing the number of macroscopic state variables, the reference empirical density 
gets approximated more accurately. We now take a closer look to the projected 
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L 


2 


3 


4 


5 


6 


7 


8 


9 


10 


P 


0.000 


0.000 


1.197- KT a 


0.840 


0.862 


0.999 


0.999 


0.999 


0.999 



Table 1. The p- values of a two-sample Kolmogorov-Smirnov test 
that compares the reference and projected empirical distributions. 
Simulation details are in the text. 



ensembles by computing the relative difference between the l-th even empirical 
moment of the projected ensemble, Ui, and the corresponding empirical moment 
of the reference ensemble U* as (Z7; — U*)/U*. Figure 1(b) shows this error as a 
function of I for different values of the number of macroscopic state variables L. 

We make two key observations. First, for I < L, the relative difference in the 
corresponding moment is of the order of the tolerance of the Newton procedure. 
This is expected, since these are the macroscopic state variables onto which the 
distribution is projected. We see that this very small error increases nevertheless 
with /; this can be explained by pointing out that the value of the moments increases 
very quickly with I, and that the equations in the Newton procedure have been 
rescaled accordingly. Second, the error in the higher moments (I > L) also decreases 
with increasing L. This indicates that the convergence of <P\l] to (p* for L — > oo is 
not only due to the fact that we project onto more moments, but also because the 
approximation of the higher order moments improves. So far, we have no complete 
theoretical justification for this observation. 

Finally, we compare the reference and projected empirical distributions using a 



two-sample Kolmogorov-Smirnov test 37 ; this classical hypothesis test results in a 



high p- value (< 1) if the two samples are likely to have been drawn from the same 
probability distribution. The results are shown in Table We clearly see a p- value 
that approaches 1 for increasing L. 



4.3.2. Error dependence on the time step. In a next experiment, we simulate up to 
time t* = 1.54, and record the microscopic state y~ at time t~ = 1.5, as well as 
the macroscopic states U\n{t~ + At) for L = 3,4,5, and the function of interest 
f p (t~ + At) for At <G [0,i* — t~\. We then project the ensemble y onto the 
macroscopic state Um(t~ + At), obtaining V[L](y~ , U^itT + At)), and denote 
the corresponding value of the stress as f p (t~ + At). We record the relative error 
\f p (t) — f p (t)\/f p (t) as a function of At. To reduce the statistical error, we report 
the averaged results of 100 realizations of this experiment. The results are shown in 
Fig. [2] where we used e = 1 in ^ . We indeed see the linear increase of the matching 
error as a function of At; notice also, as was shown in the previous experiment, that 
the matching error decreases with increasing L. From the figure, we conclude that 
the remaining statistical error is at least lower than 10 -5 . 



5. Extrapolation operator 

Next, we need to specify how the extrapolation is performed. Let the order of 
consistency be defined as follows: 

Definition 5.1 (Consistency of extrapolation). Consider a certain class of suffi- 
ciently smooth functions. An extrapolation operator £ is called consistent of order 
p e > for this class if there exist At > and C such that for all At £ [0, At ), all 
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Figure 2. Error of the stress after projecting a prior ensemble 
of FENE dumbbells onto the first L even centralized moments of 
a reference ensemble, as a function of At for several values of L. 
Simulation details are given in the text. 

n < N and all functions U in the considered class it holds 

(29) \\U n+1 - U(t n+1 )\\ < CAt p * +l , 
with 

u^^s((u(t-%; k l 00 -(At^ =0 ,5t). 

Further, we will also use the following definition of continuity: 

Definition 5.2 (Continuity of extrapolation). Consider a certain class of sufficiently 
smooth functions. An extrapolation operator £ consistent of order p e > for this 
class is called continuous for this class, if there exist Ato > and C such that for 
all At £ [0, Ato], all n < N, and all functions U\, XJi in the considered class it holds 

II* (( W' fc )) I £ ,o i ( A ^=o, st) £ (( W fc )):%„,o ; (AtOS* 5t 

(30) ^fAty n ' K 

St) 



i.fe=0,0 



In the remainder of this section, we consider two extrapolation strategies: projec- 



tive extrapolation (Subsection 5.1 ) and multistep state extrapolation (Subsection 5.2 ) 



Numerical illustrations are given in Subsection |5.3| 
5.1. Projective extrapolation. The first approach that we consider, proposed 



m 

in, k 



15 



is to extrapolate U to time t n+1 using only (some of) the time points 
t n ' k , k = 0, . . . , K , i. e., by using the sequence of points obtained in the last burst 
of microscopic simulation. If this coarse projective extrapolation is based on the 
interpolating polynomial of degree p e (with p e < K) through the parameter values 
at times t n,k , k = K — p e , . . . , K, we obtain 



(31) U n+1 = J2l s (a n )U 



K-t 
s=0 
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with 
(32) 

and the Lagrange polynomials 
(33) l s (a) = 



At, 
St 



-K, 



a(a + 1) ■ • • (a + p e ) 
s\( Pe ~sy.(-l) s (a + s)' 



Clearly, (31 1 is continuous in the sense of Definition 5.2 



Example 5.3 (Coarse projective forward Euler). The simplest, first order version 
of the above method is called coarse projective forward Euler. In this case, the 
procedure can be rewritten as 

jjn,K jjn,K—\ 



(34) 



jjn+l = jji 



,K 



(At n - K8t)U ' 



n 



st 



The procedure described above is reminiscent of a Taylor method [15] . Conse- 
quently, time integration based on this extrapolation will resemble a Taylor method 
when repeatedly extrapolating forward in time, and the global deterministic error 
will be dominated by a term of the form CAt Pe (assuming St <C At), as results from 
an accuracy analysis of coarse projective integration for deterministic microscopic 



models 57 



To assess qualitatively the influence of coarse projective extrapolation on the 
statistical error, we apply it to the linear test equation 

(35) dX{t) = aX(t) dt + bdW(t). 

Application of the one step method (f> to ( [35] ) yields 

Y n > k =R 4 {a, St, ^n.*-i)y.fc-i + s^a, b, St, rf^" 1 ), 

where and S<f, are functions depending on <f> and rj n + 1 < k - 1 are (vectors of) the 
i.i.d. random variables used by 4>. In the following, we assume that R^a^t^rf 1 ' 1 ) 
is independent of rf 1 ' 1 and can be written as R^(a,St,i] n - 1 ) = R c / ) (aSt); this holds, 
e.g., for typical Runge-Kutta methods. The above assumptions imply 



k-i 



k-i-l 



Y n > k =R 4> (aSt) k Y n ^° + J2 S*(a, b, St, rf 1 '*) R^aSt) 

i=0 

fe-1 

E Y n ' k =R 4> (a5t) k E Y n '° + E S^a, b, St, V nfi ) ^ R^aSt) 



i=0 



If we now apply the extrapolation step (31), we obtain 

Pe 

E yn+1,0 ^J2ls(a)'EY n ' K - s 



s=Q 



(36) 



K-8-1 



^2 ls(a)R<p(a6t) K ~ s J El" 1 ' + E S^a, b, St, rj n '°) ^ l s (a) ^ R^aSt) 1 



\s=0 



s=0 



i=0 



= :R E {aSt) 
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with a and l s (a) given by (32) and (33 1. In analogy to (36) we obtain also 

p„ K-s-l 

gyn+i,o = # E ( a( ft)Eyn,o + l s ( a ) Y] i? (a(5i) l E > S' <?i (a, fc, St, n n ^). 

s=0 i=0 

Thus EEy n+1 >° = Er n+1 '°, but 
VarEy ,l+1 <° = - R E (aSt) 2 VarF"' 

A"-l /min{p c ,A-l-i} 

(37) + - Var S^a, b, St, r, n >°) £ R^aStf 1 ^ i,(a) 

i=0 \ s=0 



The second summand in ( 37 1 behaves as (5ta 2pe / J for large a = At/St— K. Assuming 
At ^S> St, this results in an amplification of the statistical error with a factor At Ps /St PB 
during extrapolation. A natural question is then: how many realizations J are 
needed to obtain the same variance using a fully microscopic simulation? In that 
case, we have 

VarEy™^ +x ^jR^aSt) 2 ^^ VarF"' 

a+K-l 

(38) + Var S^, (a, MM"'°) R <t>( aSt ) 21 - 

J i=0 

Thus, for large a, the required number of realizations for a full microscopic simulation 
is smaller by a factor l/a 2p "~ 1 , i.e, J ~ n2 ^_ 1 , whereas the computational costs 
per realization increases by a factor a. This means that for large a and p e = 1 , the 
computational cost of the micro/macro acceleration technique with coarse projective 
extrapolation is similarly to that of a full microscopic simulation for a given variance. 
For large a and p e > 1, coarse projective extrapolation is even more expensive than 
a full microscopic simulation. 

To reduce statistical error, it has been proposed to use a chord based approxima- 
tion, for instance, using jj n . K - K i f or the time derivative estimate instead of JJ n,K ~ 1 



in equation (34) 46 . Instead of taking Lagrange polynomials in equation (311, we 
then have 



K 



(39) U n+1 = Y,Uan)U n < K - s , 

in which Ik{o) = 1 + K " K , Iki(oi) = ~ k-k > anc ^ ^*( a ) = otherwise (see 



Example 5.3) reduces the variance by a factor 1/(K — K\). However, the conclusion 



on the computational cost remains the same. 

5.2. Multistep state extrapolation. Because of the amplification of statistical 
error, we look into alternative extrapolation strategies. One approach, proposed 
in 53 58 , is to extrapolate U to time t n+1 using only (some of) the time points 
t %,K , i = l,...,n, i.e., by using the last point of each sequence of microscopic 
simulations, instead of a sequence of points from the last microscopic simulation. If 
this multistep state extrapolation method is based on the interpolating polynomial 
of degree p e through the parameter values at times t i,K , i = n — p e , . . . ,n, and we 
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assume equidistant coarse time steps At, we obtain 
(40) U n+1 = Y J l*{P)U n ~ s ' K , 

where 

a + K 

is the fraction of the interval At over which we extrapolate, and a and l s are defined 



as in (32 1 and (33 1. Note that such an extrapolation strategy requires a separate 
starting procedure. 

For a detailed comparison of the accuracy and stability properties of acceleration 
of the numerical integration of ODEs using projective extrapolation and multistep 



state extrapolation, we refer to 57 58 . Here, we only remark that, while the 



local error using multistep state extrapolation only differs by a factor two with 
respect to the local error of projective extrapolation, the global error is affected 
quite significantly. This is due to the fact that, when increasing At, also the time 
derivative estimate itself is taken over a larger time interval. It has been shown 
in 53 57 that this results in an error constant (see, e.g., 18 Section III. 2]) of the 
form CaAt Pe . This amplification effect will be illustrated in Subsection |5.3| 

Let us now look into the statistical error, again using the linear test equation 



(35 1. One extrapolation step of coarse multistep state extrapolation yields 



s=0 

Thus again EEY™ +1 >° = EF" +1 '°, but now 

VarEy ,i+1 <° < imaxVarY n - s ' K ^|k(0)|J • 

As fj < 1, the last factor can be bounded (independently of a). Consequently, the 
amplification of statistical error during the extrapolation does not depend on a, 
whereas the corresponding computational costs per simulation path are reduced by 
a factor a compared to a full microscopic simulation. 

5.3. Numerical illustration. We now provide a numerical result to illustrate the 
effects of extrapolation on the deterministic and statistical error. To avoid effects 



of the matching step, we consider the linear equation (27) with 02(f) = —ai(t) = 
b(t) = 1 , for which we know that macroscopic evolution closes in terms of the first 
two moments of the distribution. This microscopic SDE is discretized using an 
Euler-Maruyama scheme with St = 2 • 10~ 4 . We consider 500 realizations of a 
computational experiment with J = 1000 SDE realizations. As an initial condition, 
we sample from a standard normal distribution. We compare the sample mean 
behavior and sample standard deviation of a full microscopic simulation (which 
we will call the reference simulation) with the micro/macro acceleration algorithm 
using At = 1 ■ 10~ 3 , 2 ■ 10~ 3 , 4 • 10~ 3 , and 8 • 10~ 3 . The function of interest is chosen 
to be /(t) = E(A(t) 2 ). We denote by /(f) the approximation to the function of 
interest calculated from one realization of the reference simulation using J SDE 
realizations, and by /(f) the function of interest obtained via one realization of the 
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At = 0.008 
At = 0.004 
At = 0.002 
At = 0.001 
reference 



Figure 3. Results of micro/macro acceleration of the linear equa- 
tion (27) using L = 2 moments and projective extrapolation for 



different values of the time step At, as well as a full microscopic 
(reference) simulation. Top left: evolution of the sample means 
of the function of interest / and /. Bottom left: deterministic 
error on /. Right: evolution of the sample standard deviation of /. 
Simulation details are given in the text. 



micro/macro acceleration technique. As extrapolation techniques, we use first order 
projective extrapolation and first and second order multistep state extrapolation. 

Figure [3] shows the results for first order projective extrapolation. The left figure 
clearly shows, as expected, that the deterministic error grows with increasing At. 
However, from the right figure follows that, for an individual realization of the 
experiment, the error is dominated by the statistical error. The zoom shows that, 
for small t, the sample standard deviation grows linearly as a function of time, with 
a slope that is larger for larger At. This is in agreement with the theoretical result 
on the local propagation of statistical error. 

Next, we look at first order multistep state extrapolation, for which the results 
are shown in Fig. |4j The left figure indicates that, when comparing with projective 
extrapolation, the deterministic error grows much more rapidly with increasing At. 
On the right, we see that, while the sample standard deviation is larger than for 
the reference simulation, the sample standard deviation does not depend crucially 
on At, as is also expected from the analysis of the local propagation of statistical 
error. Note that the lower sample standard deviation for At = 8 • 10~ 3 is related to 
the fact that f(t) itself is much lower as a consequence of the large deterministic 
error (see left figure). The zoom shows that, for small t, the statistical error grows 
linearly as a function of time, with a slope that is independent of At, and is identical 
to the slope for the reference simulation. These results are in agreement with the 
theoretical result on the local propagation of statistical error. 

Finally, we consider second order multistep state extrapolation. The results are 
shown in Fig. [5] The left figure shows that the deterministic error is much better 
than for the first order version. However, the behavior of the statistical error is 
more intriguing. When zooming in to the behavior for small t, we observe that, 
over a short time interval, the sample standard deviation using the micro/macro 
acceleration technique increases at the same rate as the sample standard deviation 
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— At = 0.008 
---- At = 0.004 
At = 0.002 



Figure 4. Results of micro/macro acceleration of the linear equa- 
tion (27 1 using L = 2 moments and first order multistep state 
extrapolation for different values of the time step At, as well as a 
full microscopic (reference) simulation. Top left: evolution of the 
sample means of the stresses / and /. Bottom left: deterministic 
error on /. Right: evolution of the sample standard deviation of /. 
Simulation details are given in the text. 




V 4 




At = 0.008 
At = 0.004 
At = 0.002 
At = 0.001 
reference 



Figure 5. Results of micro/macro acceleration of the linear equa- 
tion ( 27 1 using L = 2 moments and second order multistep state 



extrapolation for different values of the time step At, as well as a 
reference fully microscopic simulation. Top left: evolution of the 
sample means of the function of interest / and /. Bottom left: 
deterministic error on /. Right: evolution of the sample standard 
deviation of /. Simulation details are given in the text. 



in the reference simulation. This corresponds to the theoretical result on the local 
propagation of statistical error. However, on longer time scales, the sample standard 
deviation for large t grows rapidly, and seems to be larger for larger At. As such, 
second order multistep state extrapolation behaves similarly to first order projective 
integration on long time scales. We suspect that the loss of this favorable error 
propagation is due to accumulation effects. This is indicated by the fact that 
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the length of the time interval on which the local theoretical results is observed 
is longer for larger values of At. Indeed, the figure indicates that the effects of 
accumulated statistical errors start to appear after a given number of extrapolations, 
independently of the size of the extrapolation step. This behavior requires additional 
analysis. 

6. Convergence results 

Using the above, we are now ready to give a definition of convergence for the 
proposed algorithm: 

Definition 6.1. Consider a sequence of restriction operators r_i 2 anc ^ 

a sequence of matching operators \P[L]) L —i 2 > an< i denote the corresponding 
numerical approximation process obtained by using 7Z[l] an d "P[V\ m Algorithm 



3.2 



by Yji], and the maximum step size by At, At = max^ =1 At n . The accelerated 
micro/macro Monte Carlo simulation is then called weakly convergent to the solution 
X of SDE ([I]) as At — > and L — > oo at any time t € I At with time order p if 
for each / e C^ p+1 \R d ,R) there exist constants At > 0, L , Cl, and Cl, with 
Cl — y for L — y oo, such that 

(42) | E f(Y [L] (t)) - E f(X(t)) \<C L + C L (Atr 

holds for all t e I At , all L > L , and all At G [0, At }. 

We first discuss convergence when extrapolation is performed as in coarse projec- 



tive integration (Subsection 6.1l. Due to the multistep nature of the extrapolation, 
proving convergence for the multistep state extrapolation method is more involved; 
Subsection 16^21 contains a result for a linear SDE. 



6.1. Convergence using projective extrapolation. The following theorem gen- 
eralizes the theorem for the convergence of one step methods due to Milstein 



(see WO 43 or also [10]). 



Theorem 6.2. Suppose the following conditions hold: 
(i) The coefficient functions a(x) and b l {x) (where b % denotes the i-th column of 
b) are continuous, satisfy a Lipschitz condition with respect to x, and belong to 
£,p+i,2(p+i)^ ^ jg,^ jg,-^ ^ _ ^ _ )TO p or non jj-Q SJJEs^ we require in addition 

that b l is differentiable and that also b l 'b l satisfies a Lipschitz condition and 

belongs to C p p 2(p+1) (I x R d , R), i=l,...,m. 
(ii) For sufficiently large r the moments E(|| Yj^j fc || 2 '') exist for k = 0, . . . ,K and 

n = 0, 1, . . . , N and are uniformly bounded with respect to L, N . 
(Hi) SDE (JT|) is uniformly weakly continuous for the set of all Y^ K . 
(iv) The one step method <p is weakly consistent of order p^,. 

(v) The sequence of matching operators is continuous for the numerical approxi- 
mation process and all sequences of macroscopic states, and consistent for all 



sequences of triples \ Y$* ,X t ' X W (t n+l >°), K [L] (Y^ 1 ' ) 



£=1,2,... 



(vi) The extrapolation is consistent of order p e > 1 and continuous for the class of 

all functions U(t) = TZ [L] {X e '" ,Y W (t)) . 
Then, the micro/macro acceleration algorithm with projective extrapolation is weakly 
convergent with time order p = \jmi{p e ,p ! f ) } . 
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Proof. Let g(s,x) := E [f(X{t n+1 )) \X{s) = xj for s G I, x G R d , and t n+1 G 7 A * 

with s < t n+1 . Due to condition (ji| g G C P ' 2 ^ >+1 ' ) 40 . Therefore, the consistency of 
</> implies that <? satisfies 

(43) |E ff (s,X t ' x (t + St)) -Eg(s,</>(t,x;5t))\ < C g {x)5t p 



3.7 



uniformly w. r.t. s £ [t°,t n+1 ] for some C g £ Cp(M d ,]R). Analogously, Corollary 
and the continuity of the matching imply that g satisfies 

|E 5 ( S ,X ti, Mi] K (t i + 1 ' ))-E ff ( S ,Yji+ 1 ' )| 

< | E g(s, X*'*' Y m (f + 1 - )) - E g L V [L] (y*£* , K [L] (X^ > Y m (f +^ )) 
+ | E g (a, V [L] (Y^,K [L] (X**^ (f^))) ) - E g (s, V [L] (Y^>°) 

< CLAt+CxWnmix^^m -K m (Y» lfi )\\. 

The last summand can be expanded as follows: 
||%](X t ' ,X '^] f (t i + 1 ' ))-% ] (l^| 1 ' )|| 

+ p [L] (X* i '°' V ii](t i + 1 ' )) -£ ('(% ] (X i< '°^l(t i ' fe )))^ o ;At i) ^ || 



rAt 



+ \\s (jn [L] (^'^V)))* ; a^, 5t) - £ ((n [L] (y$ ))* q ; At,, at 

Let §[x,](s,x;*) := ^E (#;(X(t))|X(s) = for s G /, x G R d , and f G 7 Z 

with s < t. With this definition, the continuity and consistency of the extrapolation 
imply 

||% ] (X t<,X '^(t i + 1 ' ))-% ] (Yjf 1 ' )|| 
<\\K [L] (X ti ' K ' Y ^(t i+1 '°)) - K [L] (X ti,0 ' Y m(t i+1 > )) || + C 2 (At^ +1 



+ C 3 



- ^Il^^'^^))-^,^ 

fc=l 



= || E~g [L] {?>*, Yfjf ; f + 1 ' ) - E§[L] (t*^, X t<,0 ' y W t^ 1 ' ) || + C 2 (Ai) Pe+1 



C 3 



AA P 



K 



- Eii( E .9K^ ^ (^))-E 5 K^r 

fc=i 



L 



1=1 



where we made also use of X t% '°' Y m (t i+1 '°) = ' m (* i,K )(ti+i.Q). Due to 

the consistency of <j>, altogether we obtain 

(44) 

■[£] 



lEc^X 4 *' ,Y ii] (i i+1 < )) -Eg^Y^j 1 ' )! < C L At + C , L (At) p = ( 5^ + C L ^ +1 
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uniformly w. r. t. s E [t°,t n+1 ] for some constants Cl and Cl with Cl — > for 
L — y oo . For ease of notation, in the following we will neglect the L-dependency of 
Y. Then 

Ef(X t0 - x °(t n+1 )) -Ef(Y n+1 >°) 

n K-l 

i=0 k=0 
n-1 

+ (E /(X t * ,K ^ i ' K (r+ 1 )) - E f(X ti+1 '°' Yi+1 -\t n+1 ) 



i=0 



+ Ef(X t, ' Y, (t n+1 ))-Ef(Y n+lfi ). 

Making use of X 1 ''"^'" (t n+1 ) = x t^+\x^ k '(t^ 1 ) and comb ini n g the 

last two summands we obtain 

E/(X t0 < x °(i n+1 )) -Ef(Y n+1 >°) 

= EE ( E /( xt " +1,xtI ' ' (tl ' fc+1) (t" +1 ))-E/(x tl ' fc+1 ^' fc+1 (r+ 1 ))) 

n 

+ (E /(X f '' Kyi ' K (f +1 )) - E /(X t!+1 ' ^ 1+1 '°(t" +1 ))) . 

Using X t% ' K > Y *' K (t n+1 ) = x i ' +1 '"^' Mt,5 " ,K ( t,+1 '°)(f»+ 1 ) and the definition of g this 
implies 

Ef(X t0 < x °(t n+1 )) -E/(Y n+1 <°) 

n A'-l 

= X) E (E.g(i ! '' £+1 ,X t, ' fc ^ , ' fc (^ fe + 1 )) -E« ? (^ fe+1 ,y^+ 1 

i=0 fc=0 



+ ]T (e /(X t!+1 '°^ fI ' K ^V* 1 ' )^)) - E/(X* 4+I, °' yl+1,0 (t^ 1 )) 

n A-l 

= 2 Yl (vg(t l ' k+ \X tZ ' k > Yt - k (f< k+1 )) -Eg(f' k+ \Y l < k+1 ) 

n 

+ ^(E !/ (i«+ 1 ' ,I i,,KrM! (i ,+1 ' ))^E S (f+ 1 ^r+ 1 ^ 



Thus, (|43| and (|44|) imply 

|E/(X t ^°(r+ 1 ))-E/(F [i] (t"+ 1 ))| 

n it— 1 n 

< E E VC g (Y^)5t p * +1 + J2 (C L At + C L {At) p °5t p * + C L 5t p * +1 ) , 

i=0 k=0 i=0 

which yields together with condition (JiTJ) the desired convergence. □ 
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6.2. A result using multistep state extrapolation. For multistep state extrap- 
olation, the analysis is complicated by the multistep nature of the method. In this 
subsection we therefore restrict ourselves to consider linear SDEs ( 27 1 with normally 
distributed initial values and restriction operator 



(45) W-VVarY, 

If we assume now equidistant coarse time steps and that the extrapolation step is 
given by ( 40 1 , then we obtain 



E Y n+1 -° = J2 UP) E Y n ~ s ' K , Var Y n+1 '° = ^ l s (/3) Var Y n 



s,K 
s=0 



with l s and /3 given in (J33J) and (41 1. Application of the one step method to (27 1 
yields 

Y n+1 ' K =R 4 ,{a u t n+l > K -\5t, ^+1^-1)^+1^-1 
(46) + ^(ai, aa, b, t n+l < K ~\St, rf+^ K -i) 

K-l 



11 R 4> {a 1 ,t n+1 ' k , St, 7? «+i' fc )Y n+1 '° 

k=0 

K-l K-l 

+ J2 St(a 1 ,a 2 ,b,t n+1 > k ,5t,r ] n+1 > k ) J[ R4a u t n+1 >\St, v n+1 ' 1 ), 



fc=0 i=k 

where R^ and are functions depending on <p, similar to the stability function in 
the deterministic case, and -q n+1 ' k are (vectors of) i.i.d. random variables used by 
<fi. Assuming that R^iax, t n+1 ' k , St, rj n ' k ) is independent of rf 1 ' , which holds, e.g., 
for typical Runge-Kutta methods, we obtain the multistep formulas 

K-l Pe 

EY n+l,K = -Q Jl^ ai:t n+h^ S ^J2 ls ^ EY n-s,K 

k=0 s=0 

K-l K-l 



(47) + ^^{ai,a 2 ,b,t n+l ' k ,St,rf +l ' k ) J] R^a,, t n+1 '\ St), 

k=0 i=k 
K-l p c 

YaxY n+1 ' K = Y[ R <l) (a 1 ,t n+1 > k ,St) 2 J2l s (f3)VsLrY n - s > K 

k=0 s=0 

K-l K-l 

(48) + Vtu§^a l ,a 2 ,b,t n+1 ' k ,6t,Ti n+l > k ) J] R^a,, t n+1 -\ St) 2 . 

k=0 i=k 

As due to the consistency of <fi we have R^(a±, t n+x ' , 0) = 1 and 

E ^(oi, a 2s b, t n+1 ' k ,0, r) n+1 < k ) = Var S <j> {a l ,a 2 ,b, t n+1 - k , 0, r, n+1 ' k ) = 0, 
the corresponding characteristic polynomial is given for both equations by 

P(£;/3, Pe )=^ +1 -fS s (/3)r. 

s=0 

As in the deterministic case (theory of linear multistep methods) we then obtain 
the following theorem. 
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Figure 6. Error of the stress after extrapolating and projecting 
a prior ensemble of J = 1 • 10 5 FENE dumbbells onto the first L 
even centralized moments of a reference ensemble, as a function 
of At for several values of L. Displayed is the result averaged 
over 50 realizations of the experiment. Left: First order projective 
extrapolation. Right: First order multistep state extrapolation. 
Simulation details are given in the text. 



Theorem 6.3. Assume that all roots £ of P(£; f3,p e ) — lie within the unit circle 
and that all roots with absolute value one are simple. If then the one-step method 
is weakly consistent of order p^ and given by (46 1 with R$ independent ofrj, then 
coarse multistep state extrapolation with restriction operator (45 I and extrapolation 
given by (40) is convergent of order p = min{p e ,p0} for linear SDEs (27) with 
normally distributed initial values. 



7. Numerical results 
In this section, we provide further numerical results. We first illustrate the depen- 



dence of the local error in one accelerated time step on the step size (Subsection 7.1 ). 
Subsequently, we perform a number of long-term simulations (Subsection |7.2[ ) . 

Below, we consider equation ^ in one space dimension, with F(X) the FENE 
force Q with 7 = 49, We = 1. As the velocity field, we choose n(t) = 2 • (1.1 + 
sin (7rf)), and we again sample the initial states from the invariant distribution of 
equation ^ for n(t) = 0, and use e = 1 in d9l. We discretize in time with the 
classical Euler-Maruyama scheme with time step St = 2 ■ 10~ 4 . As before, the 
macroscopic state ETt^i consists of the first L even centralized moments. We study 
the micro/macro acceleration algorithm with first order projective extrapolation, 
and with first and second order multistep state extrapolation. In all cases, we 
perform K = 1 microscopic steps before extrapolation. 

7.1. Local error. We simulate ,7=1- 10 5 realizations up to time t* — 1.6, and 
record the microscopic state y~ at time t~ = 1.4, as well as the macroscopic states 
U\L](t~ + At) for L = 3,4, 5, and the approximated function of interest f p (t~ + At) 
for At € [0,t* — t~]. We then extrapolate from time t~ to time t = t~ + At, 
and project the ensemble y onto U^{t~ + At). Subsequently, we compute the 
corresponding value of the stress as t„(t~ + At). We record the relative error 
with respect to the reference solution, \f p (t) — ? p (t)\/f p (t), as a function of At. To 
reduce the statistical error, we report the averaged results of 50 realizations of this 
experiment. Then, the statistical error has an order of magnitude of about 10~ 5 . 
The results are shown in Fig. [6] 

For projective integration, we clearly see a first order behavior as a function of At; 
this is a consequence of the amplification of the statistical error during projective 
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Figure 7. Error of the stress after extrapolating and projecting 
a prior ensemble of J = 1 • 10 3 FENE dumbbells onto the first L 
even centralized moments of a reference ensemble, as a function 
of At for several values of L. Displayed is the result averaged 
over 20 realizations of the experiment. Left: First order projective 
extrapolation. Right: First order multistep state extrapolation. 
Simulation details are given in the text. 

extrapolation. Note that, due to the presence of three competing sources of errors 
(extrapolation, matching, and statistical error), which may be of opposite signs, the 
effect of extrapolating more macroscopic state variables is not so clearly visible as 
in Fig. [2] Note that for L = 4 the statistical and matching errors seem to be a 
bit lower, such that the second order behavior of the extrapolation error is already 
apparent for the largest displayed time steps. 

For multistep state extrapolation, the situation is slightly different. Here, we 
see that, for modest gains (small At), the statistical error remains more or less 
unaffected. For At > 2 • 10~ 3 (gain factor 10), however, the error increases as 
At 2 , as a consequence of the large extrapolation error. Note that, as soon as the 
extrapolation error dominates, the error appears to be independent of the number 
of moments used. 

To emphasize the effect of the statistical error, we repeat the experiment using 
J = 1000 realizations (averaged over 20 realizations of the experiment). The results 
are shown in Fig. [7] Compared to Fig. |6j we see qualitatively the same behavior. 
Again, the error of projective extrapolation increases linearly (but it is now an 
order of magnitude larger), while, for multistep state extrapolation, larger gains 
appear to be possible since the statistical error now dominates for a wider range of 
extrapolation step sizes At. 

7.2. Long-term simulation. We now turn to a long-term simulation, and compare 
the behavior of the sample mean and sample standard deviation of a full microscopic 
simulation (which we will call the reference simulation) with the micro/macro 
acceleration algorithm. We denote by f p (t) the approximation to the function of 
interest calculated from one realization of the reference simulation using J SDE 
realizations, and by T p (t) the function of interest obtained via one realization of 
the micro/macro acceleration technique. As extrapolation techniques, we use first 
order projective extrapolation and second order multistep state extrapolation. (For 
the set-up in this example, the deterministic error of first order multistep state 
extrapolation is too high to be considered further.) 

Figure [8] (top left) shows the evolution of the stress as a function of time. We see 
that the simulation exhibits a periodic behavior, with a fast increase of the stress 
followed by a relaxation. Note that, in this problem, due to the fast variations 



ACCELERATED MICRO/MACRO MONTE CARLO SIMULATION OF SDES 



27 



1(10 




















i 















L = 2 
L = 3 
L = 4 
reference 



Figure 8. Results of micro/macro acceleration of the FENE model 
Q using At = 1 ■ and projective extrapolation for different 
numbers L of macroscopic state variables, as well as a full micro- 
scopic (reference) simulation. Top left: evolution of the sample 
means of the stresses f p and t p . Bottom left: deterministic error 
on T p . Right: evolution of the sample standard deviation of f p . 
Simulation details are given in the text. 



in n(t), there is not a very strong time-scale separation between the evolution 
of the stress tensor and the evolution of individual polymers. Hence, using the 
time-step adaptation strategy (see Subsection 3.3) will prove crucial to obtain an 
efficient algorithm. During the fast increase of the stress, we observed that the 
matching operator fails for time steps At > 1 • 1CP 3 , whereas larger accelerations are 
possible during the relaxation. Therefore, in this experiment we will use adaptive 
macroscopic time steps, as outlined in Subsection |3.3| choosing a = 0.2 and a = 1.2. 
On average, we obtained a speed-up factor of 4, meaning that microscopic simulation 
has been performed over 1/4 of the time domain. 

In a first experiment, we use At = 1 • 10~ 3 and vary the number L of macroscopic 
state variables. The results for 500 realizations of this experiment, each with an 
ensemble size of J — 5000, are shown in Fig. [8] 

We make two main observations. First, the deterministic error decreases with 
increasing L, whereas the sample standard deviation is independent of L. (The 
different lines in the plot are nearly indistinguishable.) Note also that the variance 
on the sample standard deviation is quite large in this example. Second, from Fig. [8] 
we see that the error of the micro/macro acceleration algorithm with respect to 
the reference simulation also decreases as a function of time, until it reaches a level 
of the order of the statistical error. This behavior can be attributed to the fact 
that, in this example, the macroscopic behavior of the system on long time scales 
is determined by only a few macroscopic state variables. The results for multistep 
state extrapolation (not shown) in terms of L are similar. 

In a second experiment, we fix L = 3 and consider varying At. (This experiment 
and its conclusions closely resemble the one in Subsection 5.3 ) The results for 500 
realizations, each with an ensemble size of J = 1000, are shown in Fig. [9} Figure [9] 
(bottom left) shows again that the deterministic error grows with increasing At, 
whereas the right figure illustrates that the sample standard deviation is larger 



28 



KRISTIAN DEBRABANT AND GIOVANNI SAMAEY 




Figure 9. Results of micro/macro acceleration of the FENE model 
Q using L = 3 moments and projective extrapolation for different 
values of At, as well as a full microscopic (reference) simulation. 
Top left: evolution of the sample means of the stresses f p and f p . 
Bottom left: deterministic error on f p . Right: evolution of the 



sample standard deviation of f p . 
the text. 



Simulation details are given in 
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Figure 10. Results of micro/macro acceleration of the FENE 
model ([7]) using L = 3 moments and second order multistep state 
extrapolation for different numbers L of macroscopic state variables, 
as well as a full microscopic (reference) simulation. Top left: evo- 
lution of the sample means of the stresses f p and t p . Bottom left: 
deterministic error on f p . Right: evolution of the sample standard 
deviation of f p . Simulation details are given in the text. 



for larger At. For second order multistep state extrapolation, we obtain Fig. 10 
The behavior of the sample standard deviation is similar to the linear case. When 
zooming in to the behavior for small t, we observe that, over a short time interval, the 
sample standard deviation using the micro/macro acceleration technique increases 
at the same rate as the sample standard deviation in the reference simulation, which 
corresponds to the theoretical result on the local propagation of statistical error. 
However, on longer time scales, we again see that the sample standard deviation for 
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large t grows rapidly, and seems to be larger for larger At, due to accumulation effects. 
Note that, in this case, second order multistep state extrapolation even behaves 
worse than first order projective integration on long time scales for sufficiently large 
At. 

8. Conclusions and outlook 

We presented and analyzed a micro/macro acceleration technique for the Monte 
Carlo simulation of stochastic differential equations (SDEs) in which short bursts of 
simulation using an ensemble of microscopic SDE realizations are combined with an 
extrapolation of an estimated macroscopic state forward in time. The method is 
designed for problems in which the required time step for each realization of the 
SDE is small compared to the time scales on which the function of interest evolves. 
For such systems, one often needs to take a very small microscopic time step, which 
results in a deterministic error that is much smaller than the statistical error. 

We showed that the proposed procedure converges in the absence of statistical 
error, provided the matching operator satisfies a number of natural conditions, 
and we introduced a matching operator that satisfies these conditions for Gaussian 
random variables. We also conjectured that this matching operator is suitable for 
general distributions, and provided numerical evidence to support this conjecture. 
Concerning the statistical error, a local analysis of projective extrapolation shows 
that the amplification of statistical error depends on the ratio a of macroscopic 
(extrapolation) and microscopic (simulation) time steps, while this is not the case 
for multistep state extrapolation. Numerical evidence, however, suggests that, when 
using higher order multistep state extrapolation, accumulation of statistical error 
over macroscopic time scales may nevertheless induce an a-dependent statistical 
error. 

This paper has not focused on quantifying the computational gains that can 
be expected from this method. It is clear from the description that the method 
will be more efficient when there is a bigger separation in time scales between 
the microscopic and macroscopic levels. The numerical examples in this text do 
not exhibit such a strong time-scale separation; they were mainly chosen for their 
ability to clearly illustrate the effects of the different sources of numerical error. 
However, some conclusions on efficiency can nevertheless be drawn. For a given 
required variance on the solution, the computational cost using first order projective 
extrapolation is comparable to that of a full simulation, since the former requires 
more SDE realizations due to the a-dependent amplification of statistical error. 
For first order multistep state extrapolation, extrapolation without such a drastic 
amplification of statistical error is possible, at the cost of an amplified deterministic 
error. This can be acceptable if the macroscopic function of interest changes slowly 
compared to the time step of the SDE. 

We note that, for the model problem of coupled micro/macro simulation of dilute 
polymer solutions, the amplification of statistical error using projective extrapolation 
need not be dramatic: while the computational cost of a micro/macro accelerated 
simulation and a full microscopic simulation are comparable for given variance, 
this is no longer true when coupling this Monte Carlo simulation to a PDE for the 
solvent. Indeed, when extrapolating the complete coupled system forward in time, 
a computational gain is obtained since the PDE for the solvent also does not need 
to be simulated on the whole time domain. 
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In future work, we will study stability and propagation of statistical error on long 
time scales. The numerical experiments indicate that these issues can be studied in 
a linear setting. Another open question is for which distributions the conjecture 
can be proved. From an algorithmic point of view, this work raises questions on 
the adaptive/automatic selection of all method parameters (number of moments to 
extrapolate, macroscopic time step, number of SDE realizations) to ensure a reliable 
computation with minimal computational cost. Also, a numerical comparison with 
other approaches, such as implicit approximations, could be envisaged. 
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